16S Phyloseq

Load required packages.

require(phyloseq)
 require(parallelDist)
require(ggplot2)
  require(randomForest)
require(Boruta)
  require(phyloseq)
require(stats)
  require(dplyr)
require(readr)
  require(pheatmap)
require(tidyverse)
  require(micrUBIfuns)
require(circlize)
  require(phylosmith)
require(phangorn)

Load necessary files from Dada2 processing script.

load("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/16S/DADA2outputs/dada2_outputs_16S_BeeBread_wtree.RData")

Read in metadata

metadata<-read.csv("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/BB_sample-metadata.txt", sep = "\t")
rownames(metadata)<-metadata$X
metadata$Region_Site<-paste( metadata$Region, metadata$Site, sep = "")
metadata$Treatment_Site<-paste( metadata$Treatment, metadata$Site, sep = "")

Convert files to correct format for phyloseq object and then rename samples in OTU table to match the format of the metadata file. Create phyloseq object and root phylogenetic tree for usage with distance metrics (e.g., Unifrac).

#assign object type to merge into phyloseq object
md<-sample_data(metadata)
otu<-otu_table(seqtab.nochim, taxa_are_rows = F)
tax_tab<-tax_table(taxa)
#check sample names of both otu table and metadata
#remove .fastq.gz from otu table names
sample_names(otu)<-str_split(sample_names(otu), pattern = ".fastq.gz", simplify = T) [,1]
#since they match, you can create the phyloseq object
bb16S_orig<-phyloseq(md, otu, tax_tab, fitGTR$tree)

#root tree for distance metrics like unifrac
set.seed(11)
phy_tree(bb16S_orig)<-root(phy_tree(bb16S_orig), sample(taxa_names(bb16S_orig), 1), resolve.root = TRUE)
is.rooted(phy_tree(bb16S_orig))

[1] TRUE

Pre-processing. Rename OTU IDS to something more manageable. We save the ASCV sequences in case we need them later (e.g., BLASTN searches). Create Bacterial and Archaeal phyloseq objects, but first we remove Chloroplasts.

#Extract original IDs for future reference. 
seq_df_w_OTUID<-data.frame(OTUID = 1:ntaxa(bb16S_orig), Sequence = taxa_names(bb16S_orig))
taxa_names(bb16S_orig)<-paste("OTU_" , 1:ntaxa(bb16S_orig), sep = "")

#remove chloroplast
#60 chloroplast taxa
bb16S_orig_nc<-subset_taxa(bb16S_orig, Class != "Chloroplast")

#2257 bacterial taxa - 9 non-bacteria removed. 
bb16S_bac<-subset_taxa(bb16S_orig_nc, Kingdom == "Bacteria")
#9 archaea
bb16S_arch<-subset_taxa(bb16S_orig_nc, Kingdom == "Archaea")

Remove control samples as there is only as single replicate of each.

ps_wo_control<-subset_samples(bb16S_bac, Treatment != "none")

Examine rarefaction curves. Plateau in most samples under the 5k read mark.

rarecurve(otu_table(ps_wo_control), step=50, cex=0.5)

First, check sample depths for rarefaction. Minimum sample depth is a bit over 21k reads, so we will rarefy at 20k.

bb16S_bac_rarefy

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 1957 taxa and 72 samples ]
sample_data() Sample Data:       [ 72 samples by 10 sample variables ]
tax_table()   Taxonomy Table:    [ 1957 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 1957 tips and 1956 internal nodes ]

#Exploratory analysis of Alpha diversity

Plots of Alpha Diversity by region, treatment, timepoint, and the combinations.

plot_richness(bb16S_bac_rarefy, x = "Region", measures = c("Shannon", "Observed", "Chao1"), color = "Region") + geom_boxplot() + ggtitle("Region")

plot_richness(bb16S_bac_rarefy, x = "Treatment", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Treatment")

plot_richness(bb16S_bac_rarefy, x = "Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Site") + geom_boxplot() + ggtitle("Site")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Site by timepoint")

plot_richness(bb16S_bac_rarefy, x = "Region_Treatment", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Treatment") + geom_boxplot() + ggtitle("Region_Treatment")

plot_richness(bb16S_bac_rarefy, x = "Region_Site", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Site") + geom_boxplot() + ggtitle("Region_site")

plot_richness(bb16S_bac_rarefy, x = "Region_Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Timepoint") + geom_boxplot() + ggtitle("Region_Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Treatment_Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment_Timepoint") + geom_boxplot() + ggtitle("Treatment_Timepoint")

Based upon the site and region graphs, it appears that site contains more differences than region. There are no obvious differences among treatments. The site by timepoint graph shows what appears to be an interaction of before and after treatment and site. Some sites display lower diversity following treatment and others higher, with consistent trends across all sites though.

Statistical testing of alpha diversity metrics. Using ANOVA framework and checking for assumption of normality and correlation in predictors. region and site appear to be associated, so I will drop Region for this exploratory analysis. We can revisit if we decide Region is a more accurate predictor, but it appears that from the graphs above site provides more resolution than region. https://www.pluralsight.com/guides/testing-for-relationships-between-categorical-variables-using-the-chi-square-test

#create alpha diversity table and prep data to include grouping columns
richness_BB16S_bac<-estimate_richness(bb16S_bac_rarefy, measures = c("Shannon", "Observed", "Chao1"))
#add metadata
richness_BB16S_bac$Region<-sample_data(bb16S_bac_rarefy)$Region
richness_BB16S_bac$Site<-sample_data(bb16S_bac_rarefy)$Site
richness_BB16S_bac$Timepoint<-sample_data(bb16S_bac_rarefy)$Timepoint
richness_BB16S_bac$Treatment<-sample_data(bb16S_bac_rarefy)$Treatment
#add combo columns
richness_BB16S_bac$Region_Treatment<-sample_data(bb16S_bac_rarefy)$Region_Treatment
richness_BB16S_bac$Region_Site<-sample_data(bb16S_bac_rarefy)$Region_Site
richness_BB16S_bac$Region_Timepoint<-sample_data(bb16S_bac_rarefy)$Region_Timepoint
richness_BB16S_bac$Treatment_Timepoint<-sample_data(bb16S_bac_rarefy)$Treatment_Timepoint

#test of association between 
chisq.test(richness_BB16S_bac$Region, richness_BB16S_bac$Site, simulate.p.value = TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  richness_BB16S_bac$Region and richness_BB16S_bac$Site
X-squared = 216, df = NA, p-value = 0.0004998

Tests of Shannon diversity - significant differences across sampling sites and interaction term of timepoint and site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)

    Shapiro-Wilk normality test

data:  mod$residuals
W = 0.98541, p-value = 0.5732

summary(mod)

                         Df Sum Sq Mean Sq F value   Pr(>F)    
Site                      5  6.008  1.2015  14.092 1.19e-07 ***
Treatment                 2  0.144  0.0722   0.847    0.437    
Timepoint                 1  0.130  0.1302   1.527    0.225    
Site:Treatment           10  1.383  0.1383   1.622    0.140    
Site:Timepoint            5  3.132  0.6265   7.347 7.82e-05 ***
Treatment:Timepoint       2  0.002  0.0011   0.013    0.987    
Site:Treatment:Timepoint 10  1.487  0.1487   1.743    0.108    
Residuals                36  3.070  0.0853                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Summary of model shows site to be only lower order significant predictor.
#interaction term of site x time point as well. - this was picked up in the exploratory graphs. 

#What about pairwise differences in site
TukeyHSD(mod, "Site")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$Site
               diff         lwr         upr     p adj
Ber-Bel  0.48874565  0.13009202  0.84739928 0.0028623
Dal-Bel  0.41509027  0.05643664  0.77374390 0.0154265
Sla-Bel  0.54749287  0.18883924  0.90614650 0.0006885
Spr-Bel -0.23437520 -0.59302883  0.12427843 0.3809012
Sum-Bel  0.05055529 -0.30809834  0.40920892 0.9981025
Dal-Ber -0.07365538 -0.43230901  0.28499824 0.9889939
Sla-Ber  0.05874722 -0.29990641  0.41740085 0.9961329
Spr-Ber -0.72312086 -1.08177448 -0.36446723 0.0000081
Sum-Ber -0.43819036 -0.79684399 -0.07953673 0.0092371
Sla-Dal  0.13240261 -0.22625102  0.49105623 0.8737630
Spr-Dal -0.64946547 -1.00811910 -0.29081184 0.0000529
Sum-Dal -0.36453498 -0.72318861 -0.00588135 0.0444928
Spr-Sla -0.78186808 -1.14052170 -0.42321445 0.0000018
Sum-Sla -0.49693758 -0.85559121 -0.13828396 0.0023548
Sum-Spr  0.28493049 -0.07372314  0.64358412 0.1867351

plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

#Several significant pairwise differences to accompany the significant global test for interaction term. We can
TukeyHSD(mod, "Site:Timepoint")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$`Site:Timepoint`
                    diff         lwr          upr     p adj
Ber:A-Bel:A  0.622074429  0.03364435  1.210504508 0.0305327
Dal:A-Bel:A -0.095821595 -0.68425167  0.492608484 0.9999845
Sla:A-Bel:A  0.359075661 -0.22935442  0.947505739 0.6062652
Spr:A-Bel:A -0.246621947 -0.83505203  0.341808131 0.9412092
Sum:A-Bel:A -0.192103516 -0.78053359  0.396326563 0.9903896
Bel:B-Bel:A -0.358689782 -0.94711986  0.229740297 0.6077724
Ber:B-Bel:A -0.003272907 -0.59170299  0.585157172 1.0000000
Dal:B-Bel:A  0.567312348 -0.02111773  1.155742426 0.0673126
Sla:B-Bel:A  0.377220303 -0.21120978  0.965650381 0.5353148
Spr:B-Bel:A -0.580818240 -1.16924832  0.007611838 0.0557206
Sum:B-Bel:A -0.065475688 -0.65390577  0.522954391 0.9999997
Dal:A-Ber:A -0.717896024 -1.30632610 -0.129465946 0.0067460
Sla:A-Ber:A -0.262998769 -0.85142885  0.325431310 0.9121532
Spr:A-Ber:A -0.868696377 -1.45712646 -0.280266298 0.0005130
Sum:A-Ber:A -0.814177945 -1.40260802 -0.225747867 0.0013273
Bel:B-Ber:A -0.980764211 -1.56919429 -0.392334133 0.0000701
Ber:B-Ber:A -0.625347336 -1.21377741 -0.036917258 0.0290681
Dal:B-Ber:A -0.054762081 -0.64319216  0.533667997 1.0000000
Sla:B-Ber:A -0.244854126 -0.83328420  0.343575952 0.9438927
Spr:B-Ber:A -1.202892670 -1.79132275 -0.614462591 0.0000013
Sum:B-Ber:A -0.687550117 -1.27598020 -0.099120039 0.0110389
Sla:A-Dal:A  0.454897256 -0.13353282  1.043327334 0.2675843
Spr:A-Dal:A -0.150800353 -0.73923043  0.437629726 0.9987752
Sum:A-Dal:A -0.096281921 -0.68471200  0.492148158 0.9999838
Bel:B-Dal:A -0.262868187 -0.85129827  0.325561891 0.9124155
Ber:B-Dal:A  0.092548688 -0.49588139  0.680978766 0.9999891
Dal:B-Dal:A  0.663133943  0.07470386  1.251564021 0.0162597
Sla:B-Dal:A  0.473041898 -0.11538818  1.061471976 0.2199776
Spr:B-Dal:A -0.484996645 -1.07342672  0.103433433 0.1921371
Sum:B-Dal:A  0.030345907 -0.55808417  0.618775986 1.0000000
Spr:A-Sla:A -0.605697608 -1.19412769 -0.017267530 0.0389279
Sum:A-Sla:A -0.551179177 -1.13960925  0.037250902 0.0838838
Bel:B-Sla:A -0.717765443 -1.30619552 -0.129335364 0.0067605
Ber:B-Sla:A -0.362348568 -0.95077865  0.226081511 0.5934672
Dal:B-Sla:A  0.208236687 -0.38019339  0.796666766 0.9820534
Sla:B-Sla:A  0.018144642 -0.57028544  0.606574721 1.0000000
Spr:B-Sla:A -0.939893901 -1.52832398 -0.351463823 0.0001454
Sum:B-Sla:A -0.424551349 -1.01298143  0.163878730 0.3613081
Sum:A-Spr:A  0.054518432 -0.53391165  0.642948510 1.0000000
Bel:B-Spr:A -0.112067834 -0.70049791  0.476362244 0.9999260
Ber:B-Spr:A  0.243349041 -0.34508104  0.831779119 0.9461100
Dal:B-Spr:A  0.813934295  0.22550422  1.402364374 0.0013329
Sla:B-Spr:A  0.623842250  0.03541217  1.212272329 0.0297334
Spr:B-Spr:A -0.334196293 -0.92262637  0.254233786 0.7015194
Sum:B-Spr:A  0.181146260 -0.40728382  0.769576338 0.9940289
Bel:B-Sum:A -0.166586266 -0.75501634  0.421843812 0.9970548
Ber:B-Sum:A  0.188830609 -0.39959947  0.777260687 0.9916242
Dal:B-Sum:A  0.759415864  0.17098579  1.347845942 0.0033809
Sla:B-Sum:A  0.569323819 -0.01910626  1.157753897 0.0654618
Spr:B-Sum:A -0.388714724 -0.97714480  0.199715354 0.4909324
Sum:B-Sum:A  0.126627828 -0.46180225  0.715057906 0.9997582
Ber:B-Bel:B  0.355416875 -0.23301320  0.943846953 0.6205355
Dal:B-Bel:B  0.926002130  0.33757205  1.514432208 0.0001862
Sla:B-Bel:B  0.735910085  0.14748001  1.324340163 0.0050100
Spr:B-Bel:B -0.222128458 -0.81055854  0.366301620 0.9711746
Sum:B-Bel:B  0.293214094 -0.29521598  0.881644173 0.8382101
Dal:B-Ber:B  0.570585255 -0.01784482  1.159015333 0.0643241
Sla:B-Ber:B  0.380493210 -0.20793687  0.968923288 0.5225980
Spr:B-Ber:B -0.577545333 -1.16597541  0.010884745 0.0583543
Sum:B-Ber:B -0.062202781 -0.65063286  0.526227298 0.9999998
Sla:B-Dal:B -0.190092045 -0.77852212  0.398338033 0.9911643
Spr:B-Dal:B -1.148130588 -1.73656067 -0.559700510 0.0000035
Sum:B-Dal:B -0.632788036 -1.22121811 -0.044357957 0.0259759
Spr:B-Sla:B -0.958038543 -1.54646862 -0.369608465 0.0001052
Sum:B-Sla:B -0.442695991 -1.03112607  0.145734088 0.3032108
Sum:B-Spr:B  0.515342553 -0.07308753  1.103772631 0.1334924

plot(TukeyHSD(mod, "Site:Timepoint"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon"), color = "Timepoint") + geom_boxplot() + ggtitle("Site x Timepoint")

Tests of Chao 1 - emphasizes importance of rare taxa. Site is significant as are the interactions with site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Chao1 ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)

    Shapiro-Wilk normality test

data:  mod$residuals
W = 0.9901, p-value = 0.8484

summary(mod)

                         Df Sum Sq Mean Sq F value  Pr(>F)   
Site                      5  11148  2229.6   5.007 0.00139 **
Treatment                 2    838   419.1   0.941 0.39955   
Timepoint                 1    674   674.2   1.514 0.22648   
Site:Treatment           10  11548  1154.8   2.594 0.01756 * 
Site:Timepoint            5   8414  1682.8   3.779 0.00747 **
Treatment:Timepoint       2    715   357.7   0.803 0.45570   
Site:Treatment:Timepoint 10  11692  1169.2   2.626 0.01640 * 
Residuals                36  16030   445.3                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Site and the interactions with site are significant. 

TukeyHSD(mod, "Site")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Chao1 ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$Site
              diff        lwr       upr     p adj
Ber-Bel   5.063442 -20.854212 30.981096 0.9912340
Dal-Bel -16.704018 -42.621671  9.213636 0.3961244
Sla-Bel -22.736949 -48.654603  3.180704 0.1137237
Spr-Bel -29.135890 -55.053544 -3.218236 0.0199612
Sum-Bel  -3.852226 -29.769880 22.065427 0.9975571
Dal-Ber -21.767460 -47.685114  4.150194 0.1431841
Sla-Ber -27.800392 -53.718045 -1.882738 0.0294933
Spr-Ber -34.199332 -60.116986 -8.281678 0.0041242
Sum-Ber  -8.915668 -34.833322 17.001985 0.9029200
Sla-Dal  -6.032932 -31.950585 19.884722 0.9807151
Spr-Dal -12.431872 -38.349526 13.485781 0.7011191
Sum-Dal  12.851791 -13.065862 38.769445 0.6714475
Spr-Sla  -6.398941 -32.316594 19.518713 0.9750409
Sum-Sla  18.884723  -7.032931 44.802377 0.2662758
Sum-Spr  25.283664  -0.633990 51.201317 0.0593386

plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Site x Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Site x Treatment")

Tests of Observed - unique taxa in a sample. Site is significant as are the interactions with site.

Summary of alpha diversity: Rarefying at 20k reads sufficiently captures the diversity of each sample, as demonstrated by the plateau in our rarefaction curves. Site is the most useful predictor for determining differences in alpha diversity. Several interactions with site are present for all metrics of alpha diversity. No clear trends in alpha diversity differences. Some sites display shifts responding to treatment, while others do not. Some sites show increased alpha diversity at timepoint 2 and others the opposite.

#Multivaraite exploratory analysis

Info on strata argument. Briefly, this allows us to restrict permutations to within a group and control for an effect in an analogous way to a random effect.

https://stats.stackexchange.com/questions/350462/can-you-perform-a-permanova-analysis-on-nested-data/350504#350504?newreg=8e389dc237c54d87ab240486b8b0fe33 https://stats.stackexchange.com/questions/188519/adonis-in-vegan-order-of-variables-or-use-of-strata?noredirect=1&lq=1 https://stats.stackexchange.com/questions/459407/permanova-outputs-with-or-without-random-factor https://ichthyology.usm.edu/courses/multivariate/feb_7.pdf

Both with Bray-curtis and Unifrac distances.

#extract data from phyloseq object
sd_adonis<-data.frame(sample_data(bb16S_bac_ML))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist<-phyloseq::distance(bb16S_bac_ML, method =  "bray")

#full model with no strata arguemnt. Essentiall the effect of each main effect without controlling for any potential groups. 
adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations to time points. 
adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Timepoint, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

#another way to specify the same model using adonis2. Both provide same results as would be expected. No need to run again, merely a sanity check.
#perm<-how(nperm = 10000)
#setBlocks(perm)<-with(sd_adonis, Timepoint)
#adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = perm, by = "terms")

#weighted unifrac distance
dist.uf<-phyloseq::distance(bb16S_bac_ML, method =  "wunifrac")
#full model with no strata argument, essentially the effect of each main effect without controlling for any potential groups. 
adonis2(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations to time points. 
adonis2(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Timepoint, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis2(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

Based on results, site is the main determinant of microbiome composition (r2 of 0.51 and r2 of 0.38, BC and weighted unifrac, respectively). Time point is also significant but explains much less of the variation (r2 of 0.063 and r2 of 0.021, BC and weighted unifrac, respectively). Treatment is only significant for BC (r2 of 0.024). Both remain significant for BC dissimilarity when site and timepoint are controlled for using the strata argument. However, when site is controled for using weighted unifrac dissimilarity, the only significant predictor is the site x treatment interaction. Next, I will split by time point and examine the effect of each treatment while controlling for site effects using the strata argument (strata = site).

Take away: Site is the strongest predictor of community dissimilarity. This is true for both Bray-Curtis and weighted unifrac distances.

Split by time point as these are before and after treatment.

#full models with no strata 
adonis2(dist ~  Site * Treatment, data= sd_adonis, permutations = 10000)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 10000

adonis2(formula = dist ~ Site * Treatment, data = sd_adonis, permutations = 10000)
               Df SumOfSqs      R2       F    Pr(>F)    
Site            5   6.6661 0.73148 20.9060 9.999e-05 ***
Treatment       2   0.3728 0.04091  2.9228   0.01050 *  
Site:Treatment 10   0.9264 0.10165  1.4527   0.08969 .  
Residual       18   1.1479 0.12596                      
Total          35   9.1131 1.00000                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Time point 1 (pre-treatment): Significant effect of Site (r2 of 0.73). Treatment is also a significant predictor (r2 of 0.04). This is something to note because this is supposed to be before treatment was applied. When permutations are restricted to within site, treatment and site remain significant. For weighted unifrac, treatment is insignificant (r2 of 0.025) whether site is controlled for or not.

#models with strata to control for site differences
adonis2(dist ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  strata 
Permutation: free
Number of permutations: 10000

adonis2(formula = dist ~ Site * Treatment, data = sd_adonis, permutations = 10000, strata = sd_adonis$Site)
               Df SumOfSqs      R2       F  Pr(>F)  
Site            5   5.0277 0.60488 10.5457 0.06089 .
Treatment       2   0.3306 0.03978  1.7337 0.04140 *
Site:Treatment 10   1.2372 0.14885  1.2976 0.10579  
Residual       18   1.7163 0.20649                  
Total          35   8.3119 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Time point 2 (post-treatment): Site was again a significant predictor (r2 of 0.60) and treatment was borderline significant (r2 of 0.03). However, when the effect of site was controlled for by using strata = site, treatment becomes significant (r2 of 0.039). So when controlling for site differences, the treatments contribute to differences in microbiome structure. This was not the case for weighted unifrac dissimilarity. Treatment was not a significant predictor of community dissimilarity regardless of whether or not site differences were controlled for.

Take away: Site remains the most useful predictor of community dissimilarity. Treatment effects are only observed for non-phylogentic metrics of dissimilarity when the effect of site is controlled for.

#Visualize ordinations using two methods. CAP and NMDS

#extract data from phyloseq object and run anova.cca and drop1(cca) using vegan to examine significance of constraining variables
sd_adonis<-data.frame(sample_data(bb16S_bac_ML))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
tab_adonis<-data.frame(otu_table(bb16S_bac_ML))
#run model and check global significance using anova.cca.
cca_mod<-cca(tab_adonis ~ Site + Timepoint + Treatment, sd_adonis)
anova.cca(cca_mod)

Permutation test for cca under reduced model
Permutation: free
Number of permutations: 999

Model: cca(formula = tab_adonis ~ Site + Timepoint + Treatment, data = sd_adonis)
         Df ChiSquare      F Pr(>F)    
Model     8    1.6591 3.6811  0.001 ***
Residual 63    3.5494                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#significant global model so lets look at individual parameters
drop1(cca_mod, test = "perm", permutations = 100)

          Df    AIC      F   Pr(>F)   
<none>       109.21                   
Site       5 122.00 4.6923 0.009901 **
Timepoint  1 111.59 3.9549 0.009901 **
Treatment  2 107.50 1.0161 0.455446   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

0#Site and time point are significant but treatment is not. We will not include treatment in our constrained ordination.

[1] 0

Constrained ordination using CAP

#Bray-Curtis
ord_cap<-ordinate(physeq = bb16S_bac_ML, method = "CAP", distance = "bray", formula = ~  Site + Timepoint )
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_cap, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Bray-Curtis CAP")

#strongest grouping by Site. Time point also has a smaller effect and the groupings are not nearly as clear. 

ord_cap<-ordinate(physeq = bb16S_bac_ML, method = "CAP", distance = "wunifrac", formula = ~  Site + Timepoint )
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_cap, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic()  + ggtitle("Weighted Unifrac CAP")

#strongest grouping by Site. Time point also has a smaller effect and the groupings are not nearly as clear. There are several outlier samples that ordinate far to the right on the CAP1 axis. 

Unconstrained ordination using NMDS

#Bray-Curtis
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "bray", k = 3, trymax = 1000)

Run 0 stress 0.09985802 
Run 1 stress 0.0999061 
... Procrustes: rmse 0.002736472  max resid 0.02017735 
Run 2 stress 0.1051757 
Run 3 stress 0.09985788 
... New best solution
... Procrustes: rmse 6.352421e-05  max resid 0.0003569453 
... Similar to previous best
Run 4 stress 0.1028199 
Run 5 stress 0.1020596 
Run 6 stress 0.09974139 
... New best solution
... Procrustes: rmse 0.01668137  max resid 0.09840134 
Run 7 stress 0.09985791 
... Procrustes: rmse 0.01665151  max resid 0.09813697 
Run 8 stress 0.09967763 
... New best solution
... Procrustes: rmse 0.003165557  max resid 0.02235285 
Run 9 stress 0.1020364 
Run 10 stress 0.1030262 
Run 11 stress 0.103069 
Run 12 stress 0.1020596 
Run 13 stress 0.1030691 
Run 14 stress 0.1030694 
Run 15 stress 0.1007329 
Run 16 stress 0.1055453 
Run 17 stress 0.1052402 
Run 18 stress 0.09968197 
... Procrustes: rmse 0.0008779852  max resid 0.005908349 
... Similar to previous best
Run 19 stress 0.1054156 
Run 20 stress 0.110595 
*** Solution reached

ord_nmds$stress

[1] 0.09967763

plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Bray-Curtis NMDS")

#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. 

#weighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "wunifrac")

Run 0 stress 0.07903344 
Run 1 stress 0.07903339 
... New best solution
... Procrustes: rmse 0.0003681743  max resid 0.002414965 
... Similar to previous best
Run 2 stress 0.07903322 
... New best solution
... Procrustes: rmse 0.0002776536  max resid 0.001680673 
... Similar to previous best
Run 3 stress 0.0790332 
... New best solution
... Procrustes: rmse 2.919829e-05  max resid 0.0001953481 
... Similar to previous best
Run 4 stress 0.07903329 
... Procrustes: rmse 8.114694e-05  max resid 0.0005157978 
... Similar to previous best
Run 5 stress 0.07903354 
... Procrustes: rmse 0.000149892  max resid 0.001061866 
... Similar to previous best
Run 6 stress 0.07903337 
... Procrustes: rmse 0.0002295012  max resid 0.001007456 
... Similar to previous best
Run 7 stress 0.0790333 
... Procrustes: rmse 5.741207e-05  max resid 0.0004008627 
... Similar to previous best
Run 8 stress 0.07903331 
... Procrustes: rmse 0.0001545778  max resid 0.001056759 
... Similar to previous best
Run 9 stress 0.07903354 
... Procrustes: rmse 0.000137059  max resid 0.0009598509 
... Similar to previous best
Run 10 stress 0.07903312 
... New best solution
... Procrustes: rmse 8.244484e-05  max resid 0.0005365538 
... Similar to previous best
Run 11 stress 0.07903316 
... Procrustes: rmse 7.600391e-05  max resid 0.0005325158 
... Similar to previous best
Run 12 stress 0.07903338 
... Procrustes: rmse 0.0001346872  max resid 0.0008634345 
... Similar to previous best
Run 13 stress 0.07903336 
... Procrustes: rmse 0.0001683769  max resid 0.001194867 
... Similar to previous best
Run 14 stress 0.0790334 
... Procrustes: rmse 0.0001846303  max resid 0.00131403 
... Similar to previous best
Run 15 stress 0.07903344 
... Procrustes: rmse 0.0001907853  max resid 0.00135718 
... Similar to previous best
Run 16 stress 0.0790332 
... Procrustes: rmse 0.0001058034  max resid 0.0005290044 
... Similar to previous best
Run 17 stress 0.07903351 
... Procrustes: rmse 0.0002186238  max resid 0.001555721 
... Similar to previous best
Run 18 stress 0.07903338 
... Procrustes: rmse 0.0001797126  max resid 0.001279643 
... Similar to previous best
Run 19 stress 0.07903355 
... Procrustes: rmse 0.0001992529  max resid 0.001359978 
... Similar to previous best
Run 20 stress 0.0790332 
... Procrustes: rmse 0.0001036592  max resid 0.0007355796 
... Similar to previous best
*** Solution reached

ord_nmds$stress

[1] 0.07903312

plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Weighted Unifrac NMDS")

#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. There are several outlier samples that ordinate far to the right on the NMDS1 axis. 

#Indicator taxa using Baruta Algorith Specify function for Baruta that accepts a phyloseq object.

baruta_phyloseq<-function(physeq = "Phyloseq object", GroupingVar = "Grouping variable"){
  #prep data
  df<-data.frame((otu_table(physeq)))
  md<-data.frame(sample_data(physeq))
  tax<-data.frame(tax_table(physeq))
  rownames(tax)<-str_replace(rownames(tax), pattern = "[[.]]", replacement = "=")
  #run simper
 
                  boruta_layer<-Boruta(df , factor(md[,GroupingVar]), 
                      doTrace = 1, ntree = 1000, maxRuns = 1000) # adjust factor!!!!

                  B1<-data.frame(boruta_layer$finalDecision)  
                  B2<-data.frame(apply(boruta_layer$ImpHistory, 2, mean))
                  names(B1)<-"Layer_decision"
                  names(B2)<-"Layer_importance"
                  BB1<-merge(B1, B2, by="row.names", all.x=TRUE)
                  #sorting the results
                  BB1_sort<-BB1[order(BB1$Layer_importance, decreasing = TRUE),] ## this is to give an overview and choose cutoff
                  #picking up important with cutoff
                  BB1_sort_cutoff<-filter(BB1, Layer_importance >= 2) ## number is your importance cutoff

                        #picking up OTUs from full otutable
                        signOTU<-as.matrix(BB1_sort_cutoff$Row.names) #vector with OTUs
                        signOTU<-str_replace(signOTU, pattern = "[[.]]", replacement = ".")
                        #my_taxa
                        tax$OTU <-rownames(tax)
                        signOTU<-str_replace(signOTU, pattern = "[[.]]", replacement = "=")
                        signmy_taxa<-filter(tax, OTU %in% signOTU) 
                        signmy_taxa<-as.data.frame(signmy_taxa)
                        taxa_keep<-signmy_taxa$OTU
return(taxa_keep)
}

Run Baruta algorith with full data set. We are looking for taxa that explain differences among sites.

baruta_Site_indicators<-baruta_phyloseq(physeq = bb16S_bac_ML, GroupingVar = "Site")
baruta_Site_indicators_ps_object<-subset_taxa(bb16S_bac_ML, taxa_names(bb16S_bac_ML) %in% baruta_Site_indicators)

https://jokergoo.github.io/ComplexHeatmap-reference/book/ Heatmap to visualize differnces with all samples together.

#create subset phyloseq object and pull out taxonomic info
baruta_Site_indicators_ps_object_tax<-tax_table(baruta_Site_indicators_ps_object)
baruta_Site_indicators_ps_object_tax<-data.frame(baruta_Site_indicators_ps_object_tax)
#rename NA as unknown
baruta_Site_indicators_ps_object_tax$Class[is.na(baruta_Site_indicators_ps_object_tax$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax$Family[is.na(baruta_Site_indicators_ps_object_tax$Family)]<-"Unknown"
baruta_Site_indicators_ps_object@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object)[,3]))

#create color pallete. This goes from green to white depending on z.score of individual taxa
col_fun<-colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
#Heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with otuIDs
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object), row_dend_reorder =TRUE, col= col_fun)

Based off of initial figure, it looks like there is a strong grouping based on Baruta ID’ed taxa, and there is pattern for every other row. This indicates differences in time points within each sites. We might want to split by time point and then repeat analysis.

Since we observe what appears to be a time point signature, we can split by time point and rerun the Baruta algorithm.

#pre-treatment samples
bb16S_bac_ML_A<-subset_samples(bb16S_bac_ML, Timepoint == "A")
#run algo and subset phyloseq object
baruta_Site_indicators_A<-baruta_phyloseq(physeq = bb16S_bac_ML_A, GroupingVar = "Site")
baruta_Site_indicators_ps_object_A<-subset_taxa(bb16S_bac_ML_A, taxa_names(bb16S_bac_ML_A) %in% baruta_Site_indicators_A)
# pull out taxonomic info
baruta_Site_indicators_ps_object_tax_A<-tax_table(baruta_Site_indicators_ps_object_A)
baruta_Site_indicators_ps_object_tax_A<-data.frame(baruta_Site_indicators_ps_object_tax_A)
#rename NA as unknown
baruta_Site_indicators_ps_object_tax_A$Class[is.na(baruta_Site_indicators_ps_object_tax_A$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_A$Family[is.na(baruta_Site_indicators_ps_object_tax_A$Family)]<-"Unknown"
#pull out taxonomy table for naming purposes
baruta_Site_indicators_ps_object_A@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_A)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_A)[,3]))

#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_A, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_A)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with otuIDs
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_A, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_A), row_dend_reorder =TRUE, col= col_fun)

Strong grouping of sites based on Baruta ID’ed taxa at the first sampling.

#post-treatment samples
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#run algo and subset phyloseq object
baruta_Site_indicators_B<-baruta_phyloseq(physeq = bb16S_bac_ML_B, GroupingVar = "Site")
baruta_Site_indicators_ps_object_B<-subset_taxa(bb16S_bac_ML_B, taxa_names(bb16S_bac_ML_B) %in% baruta_Site_indicators_B)
#extract taxonomy info and rename NAs as unknown
baruta_Site_indicators_ps_object_tax_B<-tax_table(baruta_Site_indicators_ps_object_B)
baruta_Site_indicators_ps_object_tax_B<-data.frame(baruta_Site_indicators_ps_object_tax_B)
baruta_Site_indicators_ps_object_tax_B$Class[is.na(baruta_Site_indicators_ps_object_tax_B$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B$Family[is.na(baruta_Site_indicators_ps_object_tax_B$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B)[,3]))

#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Again, strong groupings of sites by Baruta ID’ed taxa at time point 2.

Next, we will look for taxa that consistently explain the differences among treatments following application. Again we will use the Baruta algorithm, but this time, we will specify the grouping variable as Treatment.

Merge by treatments to see if there are consistencies across treatments. No significant indicators.

#subset to post-treatment samples only
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#summarize by Treatment_Site
bb16S_bac_ML_B_merged<-merge_samples(bb16S_bac_ML_B, "Treatment")
sample_data(bb16S_bac_ML_B_merged)$Treatment<-rownames(sample_data(bb16S_bac_ML_B_merged))
#run algo and filter phyloseq object
bb16S_bac_ML_B_merged<-prune_taxa(taxa_sums(bb16S_bac_ML_B_merged) > 0, bb16S_bac_ML_B_merged)
baruta_Site_indicators_B_merged<-baruta_phyloseq(physeq = bb16S_bac_ML_B_merged, GroupingVar = "Treatment")
#no taxa indicative of treatment x site interaction.
#Stop this inquiry here. Do not run below because no taxa are identified. 

baruta_Site_indicators_ps_object_B_merged<-subset_taxa(bb16S_bac_ML_B_merged, taxa_names(bb16S_bac_ML_B_merged) %in% baruta_Site_indicators_B_merged)
#extract taxonomy and rename NA to unknown. 
baruta_Site_indicators_ps_object_tax_B_merged<-tax_table(baruta_Site_indicators_ps_object_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged<-data.frame(baruta_Site_indicators_ps_object_tax_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged$Class[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B_merged$Family[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B_merged@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B_merged)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B_merged)[,3]))
#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Merge by treatment site and rerun algo. No significant indicators.

#subset to post-treatment samples only
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#summarize by Treatment_Site
bb16S_bac_ML_B_merged<-merge_samples(bb16S_bac_ML_B, "Treatment_Site")
sample_data(bb16S_bac_ML_B_merged)$Treatment_Site<-rownames(sample_data(bb16S_bac_ML_B_merged))
#run algo and filter phyloseq object
bb16S_bac_ML_B_merged<-prune_taxa(taxa_sums(bb16S_bac_ML_B_merged) > 0, bb16S_bac_ML_B_merged)
baruta_Site_indicators_B_merged<-baruta_phyloseq(physeq = bb16S_bac_ML_B_merged, GroupingVar = "Treatment_Site")
#no taxa indicative of treatment x site interaction.
#Stop this inquiry here. Do not run below becasue no taxa are identified. 
baruta_Site_indicators_ps_object_B_merged<-subset_taxa(bb16S_bac_ML_B_merged, taxa_names(bb16S_bac_ML_B_merged) %in% baruta_Site_indicators_B_merged)
#extract taxonomy and rename NA to unknown. 
baruta_Site_indicators_ps_object_tax_B_merged<-tax_table(baruta_Site_indicators_ps_object_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged<-data.frame(baruta_Site_indicators_ps_object_tax_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged$Class[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B_merged$Family[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B_merged@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B_merged)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B_merged)[,3]))
#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Exploratory tree of indicators from initial Baruta and then split by timpeoints.

plot_tree(baruta_Site_indicators_ps_object, color="Site", shape = "Timepoint", label.tips = "taxa_names")

plot_tree(baruta_Site_indicators_ps_object_A, color="Site", shape = "Treatment", label.tips = "taxa_names")

plot_tree(baruta_Site_indicators_ps_object_B, color="Site", shape = "Treatment", label.tips = "taxa_names")

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

#Core and unique taxa

Core take aways: The core was defined as those taxa occuring in atleast 90% of samples and having a relative abudnance of greater than 0.00001. When collapsed at the genus level, there are 33 members in the core. Using ASVs, the core consists of 49 members. Looking at the phylogenetic trees, the members of Rickettsiales condense into a single tip at the genus level. Blast-ing this sequence returns plant DNA. I wonder if this is a common plant visited by these bees? Either way, the core accoutns for a mean 79% +/- 0.007 (sd) of the total reads in each sample.

Taxonomic composition of core: Phylum: 20 proteobacteria, 8 firmicutes, 1 cyanobacteria, 2 bacteriodota, 2 actinobacteira.

Class: 12 members of gammaproteobacteria, 8 of alphaproteobacteria, 6 bacilli, 2 bacteriodia, clostridia, and actinobacteria, 1 cynanobacteria.

#Determine those taxa unique to a site at Genus level
unique_genera_to_region_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tax_genus, treatment = "Region")
unique_genera_to_site_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tax_genus, treatment = "Site")
#asv level
unique_asv_to_region_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy, treatment = "Region")
unique_asv_to_site_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy, treatment = "Site")

#genera
#186 genera are unique to one region. 176 genera are unique to a site. This means that only 10 Taxa unique to a region are not shared for the sites within that region. 
unique_genera_to_site_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_site_tax))
unique_genera_to_region_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_region_tax))

#asv
#1510 asvs are unique to one region. 1484 asvs are unique to a site. This means that only 26 asvs unique to a region are not shared for the sites within that region. 
unique_asv_to_site_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_asv_to_site_tax))
unique_asv_to_region_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_asv_to_region_tax))

plot_bar(unique_genera_to_site_tax_ps, x = "Region_Site", fill = "Genus") + scale_fill_manual(values = col_vector) +geom_bar(stat="identity")

plot_bar(unique_genera_to_region_tax_ps, x = "Region", fill = "Genus") + scale_fill_manual(values = col_vector) +geom_bar(stat="identity")

#plot treees
plot_tree(unique_genera_to_site_tax_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )

plot_tree(unique_genera_to_region_tax_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )

#percentage of reads accounted for by core. 
summary(sample_sums(unique_genera_to_site_tax_ps)/25000)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000000 0.0000800 0.0002600 0.0006411 0.0007300 0.0086400 

summary(sample_sums(unique_genera_to_region_tax_ps)/25000)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000000 0.0000800 0.0002800 0.0006817 0.0007700 0.0086400 

#taxonomic composition
table(tax_table(core_genus)[,2])

Actinobacteriota     Bacteroidota    Cyanobacteria       Firmicutes   Proteobacteria 
               2                2                1                8               20 

table(tax_table(core_genus)[,3])

     Actinobacteria Alphaproteobacteria             Bacilli         Bacteroidia          Clostridia      Cyanobacteriia Gammaproteobacteria 
                  2                   8                   6                   2                   2                   1                  12 

#pie charts of endemic taxa taxonomy
#Filter to regions and create pie charts for each. I've chosen to plot regions becasue its a bit simplier but we could quickly change this to sites. 
#CP
CP<-subset_samples(unique_genera_to_site_tax_ps, Region == "CP")
CP<-filter_taxa(CP, function(x) sum(x) >= 1 , TRUE)
pie_data_CP<-data.frame(table(tax_table(CP)[,2]))
names(pie_data_CP) <- c("Class", "Frequency")
ggplot(pie_data_CP, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("CP endemic genera at Phylum")

#NE
NE<-subset_samples(unique_genera_to_site_tax_ps, Region == "NE")
NE<-filter_taxa(NE, function(x) sum(x) >= 1 , TRUE)
pie_data_NE<-data.frame(table(tax_table(NE)[,2]))
names(pie_data_NE) <- c("Class", "Frequency")
ggplot(pie_data_NE, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("NE endemic genera at Phylum")

#SE
SE<-subset_samples(unique_genera_to_site_tax_ps, Region == "SE")
SE<-filter_taxa(SE, function(x) sum(x) >= 1 , TRUE)
pie_data_SE<-data.frame(table(tax_table(SE)[,2]))
names(pie_data_SE) <- c("Class", "Frequency")
ggplot(pie_data_SE, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("SE endemic genera at Phylum")

#WV
WV<-subset_samples(unique_genera_to_site_tax_ps, Region == "WV")
WV<-filter_taxa(WV, function(x) sum(x) >= 1 , TRUE)
pie_data_WV<-data.frame(table(tax_table(WV)[,2]))
names(pie_data_WV) <- c("Class", "Frequency")
ggplot(pie_data_WV, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("WV endemic genera at Phylum")

Endemic taxa take aways: 186 genera are unique to one region. 176 genera are unique to a site. This means that only 10 Taxa unique to a region are not shared for the sites within that region. Endemic taxa are currently defined as those occurring only at a single site. This creates its own set of issues in that many taxa are only found in a single sample. I’ve reached out to the PhyloSmith package creator to ask about augmenting this funciton to relax this a bit to include those taxa occurring in two or three sites. Endemic taxa are much less abundant and account for a max of 0.0086 % of the reads at a site (mean 0.00064%).

Alternate approach to core and endemic taxa based on phylognetic groupings as opposed to collpasing by taxonomy. Shared and unique taxa among sites by phylogeny

#cut tree at two differnt levels for exploratory purposes.
bb16S_bac_rarefy_tip_merged0_1<-tip_glom(physeq = bb16S_bac_rarefy, h = 0.1)
plot_tree(bb16S_bac_rarefy_tip_merged0_1, label.tips = "Family", text.size = 3, color = "Region")
#275 members
bb16S_bac_rarefy_tip_merged0_2<-tip_glom(physeq = bb16S_bac_rarefy, h = 0.2)
plot_tree(bb16S_bac_rarefy_tip_merged0_2, label.tips = "Family", text.size = 3, color = "Region")
#core - 24 members in core 
core_genus_tip<-taxa_core(phyloseq_obj = bb16S_bac_rarefy_tip_merged0_2, treatment = "Site", frequency = 0.90, abundance_threshold = 0.00001 ) 
#plot tip core
plot_tree(core_genus_tip, color = "Site", shape = "Timepoint", label.tips = "Family", text.size = 3 )


#Determine those taxa unique to a site in 0.2 cut tree
unique_genera_to_region_tip<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tip_merged0_2, treatment = "Region")
unique_genera_to_site_tip<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tip_merged0_2, treatment = "Site")

# 137 unique to a site and 144 unique to a region. 
unique_genera_to_site_tip_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_site_tip))
plot_tree(unique_genera_to_site_tip_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )
unique_genera_to_region_tip_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_region_tip))
plot_tree(unique_genera_to_region_tip_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )

Relevant bee bread microbiome literature:

DOI: https://doi.org/10.1128/genomeA.00247-18